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ABSTRACT 

We present FitSKIRT, a method to efficiently fit radiative transfer models to UV/optical images of dusty galaxies. These images 
have the advantage that they have better spatial resolution compared to FIR/submm data. FitSKIRT uses the GAlib genetic algorithm 
library to optimize the output of the SKIRT Monte Carlo radiative transfer code. Genetic algorithms prove to be a valuable tool in 
handling the multi- dimensional search space as well as the noise induced by the random nature of the Monte Carlo radiative transfer 
code. FitSKIRT is tested on artificial images of a simulated edge-on spiral galaxy, where we gradually increase the number of fitted 
parameters. We find that we can recover all model parameters, even if all 1 1 model parameters are left unconstrained. Finally, we apply 
the FitSKIRT code to a V-band image of the edge-on spiral galaxy NGC4013. This galaxy has been modeled previously by other 
authors using different combinations of radiative transfer codes and optimization methods. Given the different models and techniques 
and the complexity and degeneracies in the parameter space, we find reasonable agreement between the different models. We conclude 
that the FitSKIRT method allows comparison between different models and geometries in a quantitative manner and minimizes the 
need of human intervention and biasing. The high level of automation makes it an ideal tool to use on larger sets of observed data. 
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1. Introduction 

In the past few years, the importance of the interstellar dust 
medium in galaxies has been widely recognized. Dust regulates 
the physics and chemistry of the interstellar medium, and re- 
processes about one third of all stellar emission to infrared and 
submm emission. Nevertheless, the amount, spatial distribution 
and physical properties of the dust grains in galaxies are hard to 
nail down. 

The most direct method to trace the dust grains in galax- 
ies is to measure the thermal emission of the dust grains at 
mid-infrared (MIR), far-infrared (FIR) and submm wavelengths. 
Obtaining total dust masses from MIR, FIR or submm images 
should in principle be straightforward, as dust emission is typ- 
ically optically thin at these wavelengths and thus total dust 
masses can directly be estimated from the observed spectral en- 
ergy distribution. There are a number of problems with this ap- 
proach, however. One complication is the notoriously uncertain 
value of the dust emissivity at long wavelengths (Bianchi et al. 
2003; Kramer et al. 2003; Alton et al. 2004; Shirley et al. 2011) 
and the mysterious anti-correlation between the dust emissiv- 
ity index and dust temperature (Dupac et al. 2003; Shetty et al. 
2009; Veneziani et al. 2010; Juvela & Ysard 2012b,a; Kelly 
et al. 2012; Smith et al. 2012a). A second problem is the dif- 
ficulty to observe in the FIR/submm window, which necessar- 
ily needs to be done from space using cryogenically cooled in- 
struments. Until recently, the available FIR instrumentation was 
characterized by limited sensitivity and wavelength coverage, 
and the submm window was largely unexplored. This has now 
partly changed thanks to the launch of the Herschel (Pilbratt 
et al. 2010) satellite, but also this mission has a very limited 
lifetime. Finally, the third and most crucial limitation is the poor 
spatial resolution of the available FIR/submm instruments (typ- 
ically tens of arcsec). This drawback is particularly important 



if we want to determine the detailed distribution of the dust in 
galaxies rather than just total dust masses. This poor spatial res- 
olution limits a detailed study of the dust medium to the most 
nearby galaxies only (e.g. Meixner et al. 2010; Bendo et al. 
2010, 2012; Xilouris et al. 2012; Foyle et al. 2012; De Looze 
et al. 2012b; Fritz et al. 2012; Smith et al. 2012a; Galametz et al. 
2012). Moreover, several authors have demonstrated that even 
estimating total dust masses from global fluxes induces an error 
due to resolution effects (Galliano et al. 2011; Galametz et al. 
2012). 

The alternative method to determine the dust content in 
galaxies is to use the extinction effects of the dust grains on the 
stellar emission in the UV, optical or near-infrared (NIR) win- 
dow. This wavelength range has the obvious advantages that ob- 
servations are very easy and cheap, and that the spatial resolution 
is an order of magnitude better than in the FIR/submm window. 
Furthermore, the optical properties of the dust are much better 
determined in the optical than at FIR/submm wavelengths. The 
main problem in using this approach is the difficulty to trans- 
late attenuation measurements from broadband colors to actual 
dust masses, mainly because of the complex and often counter- 
intuitive effects of the star/dust geometry and multiple scatter- 
ing (Disney et al. 1989; Witt et al. 1992; Byun et al. 1994; di 
Bartolomeo et al. 1995; Baes & Dejonghe 2001b; Inoue 2005). 
Simple recipes that directly link an attenuation to a dust mass 
are clearly not sufficient; the only way to proceed is to per- 
form detailed radiative transfer calculations that do take into 
account the necessary physical ingredients (absorption, multi- 
ple anisotropic scattering) and that can accommodate realistic 
geometries. Fortunately, several powerful 3D radiative transfer 
codes with these characteristics have recently been developed 
(e.g. Gordon et al. 2001; Baes et al. 2003, 2011; Jonsson 2006; 
Bianchi 2008; Robitaille 2011). 
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A state-of-the-art radiative transfer code on itself, however, 
is not sufficient to determine the dust content of galaxies from 
UV/optical/NIR images. A radiative transfer simulation typi- 
cally starts from a 3D distribution of the stars and the dust in 
a galaxy model and calculates how this particular system would 
look for an external observer at an arbitrary viewing point, i.e. 
it simulates the observations. The inverted problem, determin- 
ing the 3D distribution of stars and dust from a given reference 
frame, is a much harder nut to crack. It requires the combination 
of a radiative transfer code and an optimization procedure to con- 
strain the input parameters. Several fitting codes that combine a 
radiative transfer code with an optimization algorithm have been 
set up (Xilouris et al. 1997; Steinacker et al. 2005; Robitaille 
et al. 2007; Bianchi 2007; Schechtman-Rook et al. 2012). All 
too often, however, this optimization procedure is neglected and 
chi-by-eye models are presented as reasonable alternatives. 

To optimize this given problem it is important to realize 
that the parameter space is quite large, easily going up to 10 
free parameters or more. As discussed before, the complexity of 
the dust/star geometry and scattering off dust particles is often 
counter-intuitive and results in a non-linear, non-differentiable 
search space with multiple local extrema . Furthermore, there is 
another inconvenient feature induced by the radiative transfer 
code. Most state-of-the-art radiative transfer codes are Monte 
Carlo codes where the images are constructed by detecting a 
number of predefined photon packages. Because of the intrin- 
sic randomness of the Monte Carlo code, images always contain 
a certain level of Poisson noise. If one runs a forward Monte 
Carlo radiative simulation, it is usually straightforward to make 
the number of photon packages in the simulation so large that 
this noise becomes negligibly small. However, if we want to cou- 
ple the radiative transfer simulation to an optimization routine, 
typically many thousands of individual simulations need to be 
performed, which inhibits excessively long run times for each 
simulation and hence implies that there will always be some 
noise present. As a result, we have to fit noisy models to noisy 
data, as reduced CCD data will always contain a certain level of 
noise (Newberry 1991). 

For complex optimization problems like this, it is not rec- 
ommended to use classical optimization methods like the down- 
hill simplex or Levenberg-Marquardt methods, but rather apply 
a stochastic optimization method instead. The main advantage of 
stochastic methods such as simulated annealing, random search, 
neural networks, and genetic algorithms over most classical op- 
timization methods is their ability to leave local extrema and 
search over broad parameter spaces (Fletcher et al. 2003; Theis 
& Kohle 2001; Rajpaul 2012b). Genetic algorithms (Goldberg 
1989; Mitchell 1998) are one class of the stochastic optimization 
methods that stands out when it comes to noisy handling because 
it works on a set of solutions rather than iteratively progress- 
ing from one point to another. Genetic algorithms are often ap- 
plied to optimize noisy objective functions (Metcalfe et al. 2000; 
Larsen & Humphreys 2003). During the past decade, genetic al- 
gorithms have become increasingly more popular in numerous 
applications in astronomy and astrophysics ranging from cos- 
mology and gravitational lens modeling to stellar structure and 
spectral fitting (Charbonneau 1995; Metcalfe et al. 2000; Theis 
& Kohle 2001; Larsen & Humphreys 2003; Fletcher et al. 2003; 
Liesenborgs et al. 2006; Baier et al. 2010; Schechtman-Rook 
et al. 2012; Rajpaul 2012a). For a recent overview of the use of 
genetic algorithms in astronomy and astrophysics, see Rajpaul 
(2012b). 

In this paper, we present FitSKIRT, a tool designed to effi- 
ciently model the stellar and dust distribution in dusty galaxies 



by fitting radiative transfer models to UV/optical/NIR images. 
The optimization routine behind the FitSKIRT code is based 
on genetic algorithms, which implies that the code is able to 
efficiently explore large, complex parameter spaces and conve- 
niently handle the noise induced by the radiative transfer code. 
Bias by human intervention is kept to an absolute minimum and 
the results are determined in an objective manner. In Section 2 
the radiative transfer code and the genetic algorithm are dis- 
cussed in more detail. We also apply the genetic algorithms li- 
brary on a complex but analytically tractable optimization prob- 
lem to check its reliability and efficiency. This knowledge is then 
used to explain the major steps and features of the FitSKIRT pro- 
gram. In Section 3 we test the FitSKIRT program on reference 
images of which the input values are exactly known. This should 
allow us to have a closer look at the complexity of the problem 
and FitSKIRTs ability to obtain reasonable solutions in an ob- 
jective way. In Section 4 we apply FitSKIRT to determine the 
intrinsic distribution of stars and dust in the galaxy NGC4013, 
and we compare the resulting model to similar models obtained 
by Bianchi (2007) and Xilouris et al. (1999). Section 5 sums up. 

2. FitSKIRT 

2.1. The SKIRT radiative transfer code 

SKIRT (Baes et al. 2011) is a 3D continuum Monte Carlo radia- 
tive transfer code, initially developed to investigate the effects 
of dust absorption and scattering on the observed gas and stellar 
kinematics of dusty galaxies (Baes & Dejonghe 2002; Baes et al. 
2003). The code has continuously been adapted and upgraded to 
a general and multi-purpose dust radiative transfer code. It now 
includes many advanced techniques to increase the efficiency, in- 
cluding forced scattering (Mattila 1970; Witt 1977), the peeling- 
off technique (Yusef-Zadeh et al. 1984), continuous absorption 
(Lucy 1999), smart detectors (Baes 2008) and frequency distri- 
bution adjustment (Bjorkman & Wood 2001; Baes et al. 2005). 
The code is capable of producing simulated images, spectral en- 
ergy distributions, temperature maps and observed kinematics 
for arbitrary 3D dusty systems. SKIRT is completely written in 
C++ in an object-oriented programming fashion so it is effort- 
less and straightforward to implement and use different stellar or 
dust geometries, dust mixtures, dust grids, etc. 

The default mode in which SKIRT operates is the panchro- 
matic mode. In this mode, the simulation covers the entire wave- 
length regime from UV to millimeter wavelengths and guar- 
antees an energy balance in the dust, i.e. at every location in 
the system, the emission spectrum of the dust at infrared and 
submm wavelengths is determined self-consistently from the 
amount of absorbed radiation at UV/optical/NIR wavelengths. 
In this mode, the SKIRT code is easily parallelized to run on 
shared memory machines using the OpenMP library: every sin- 
gle thread deals with a different wavelength. This panchromatic 
SKIRT code has been used to predict and interpret the far- 
infrared emission from a variety of objects, including edge-on 
spiral galaxies (Baes et al. 2010, 2011; De Looze et al. 2012b), 
elliptical galaxies (De Looze et al. 2010; Gomez et al. 2010), ac- 
tive galactic nuclei (Stalevski et al. 2011, 2012) and post-AGB 
stars (Vidal & Baes 2007). 

Besides the panchromatic mode, it is also possible to run 
SKIRT in a monochromatic mode so it is less time consum- 
ing when one wants to model images at one particular wave- 
length. In this mode, the OpenMP parallelization is included by 
distributing the different photon packages in the simulation over 
all available threads, which implies that some precautions must 
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Fig. 1. The function defined by equation (1) for n = 9. It contains many steep local maxima surrounding the global maximum at 
(x, y) = (0.5, 0.5) and is therefore a strong test function for global optimization algorithms. 



be taken to avoid race conditions. This monochromatic mode is 
particularly useful when one wants to fit radiative transfer mod- 
els to a particular observed image at UV, optical or near-infrared 
wavelengths, which we deal with in this paper. 

2.2. The GAlib genetic algorithm library 

Genetic algorithms are problem solving systems based on evo- 
lutionary principles. In essence evolution theory describes an 
optimization process of a population to a given environment. 
The core difference between genetic algorithms and most other 
stochastic methods is that a genetic algorithm works with a set 
(population) of possible solutions (individuals) to the problem 
(environment). Each individual consists of a number of param- 
eters (genes). For each of these genes there is a number of pos- 
sible values which we call alleles. These alleles do not have to 
be a discrete set, but can be defined as a range or pool where the 
genes should be drawn from. 

The algorithm starts by defining both the size and the con- 
tent of the initial population (generation 0). The individuals can 
be created randomly from the gene pool or they can be manu- 
ally defined. Each of the initial solutions is then evaluated and 
given a certain "fitness" value. The individuals that meet cer- 
tain criteria of fitness are then used to crossover and produce 
the first offspring (generation 1). Another, more convenient way 
of determining which individuals are fit for reproduction is by 
determining a crossover rate. If, for example, the crossover rate 
is set to 0.6, the 60% best individuals will be used to produce 
offspring. We can also define a mutation rate in a similar way. 
This rate determines what the odds are for a certain gene (and 
not individuals) of the offspring to undergo a random mutation. 
In practice this is done by removing one of the gene values and 
replacing it with a new possible value. After this step we return 
the fitness value. The individuals fit for reproduction are then se- 
lected and a new generation (generation 2) is created. This cycle 
is then repeated until a certain fitness value is obtained or until 
a pre-defined number of generations is reached. Since the bet- 
ter individuals of a population are always preferred, we expect 
the population average to shift to a better value where the better 
individuals are preferred again, etc. 



Genetic algorithms have been applied successfully to a large 
range of test problems, and are nowadays widely used as a reli- 
able class of global optimizers. They are becoming increasingly 
popular as a tool in various astrophysical applications. Most 
astrophysical applications so far have used the publicly avail- 
able PIKAIA code (Charbonneau 1995), originally developed in 
Fortran 77 and now available in Fortran 90 and IDL as well. 

We have selected the publicly available GAlib library (Wall 
1996) for our purposes. One of the main reasons for choos- 
ing GAlib is that, like the SKIRT radiative transfer code, it is 
written in C++ and thus guarantees straightforward interfacing. 
Moreover, the library has been properly tested and adapted over 
the years and has been applied in many scientific and engineer- 
ing applications. It comes with an extensive overview of how 
to implement a genetic algorithm and examples illustrating cus- 
tomizations of the GAlib classes. To illustrate the performance 
of the GAlib library, we have applied it to the same test function 
as first suggested by Charbonneau (1995) but also investigated 
by others (Canto et al. 2009; Rajpaul 2012b). The goal of the 
exercise is the find the global maximum of the function 

fix, y) = [ 1 6 x ( 1 - x) y ( 1 - y) sininnx) sin(^7ry) j ( 1 ) 

in the unit square with x and y between and 1 . The function 
is plotted in Figure 1 for n = 9. It is clear that this function is 
a severe challenge for most optimization techniques: the search 
space contains many steep local maxima surrounding the global 
maximum at (x,y) = (0.5,0.5). In addition, the global maxi- 
mum is not significantly higher compared to the surrounding lo- 
cal maxima. A classical hill climbing method would most defi- 
nitely get stuck in a local maximum. An iterative hill climbing 
method, which restarts at randomly chosen points, would also 
be able to solve this problem but it would normally take much 
longer to do so compared to a genetic algorithm. 

We demonstrate how the GAlib genetic algorithm library be- 
haves in this complex search space by plotting the contours and 
overplotting the position of the individuals for each generation. 
The results are shown in Figure 2. For these tests we use the 
same values for the mutation probability and crossover rate as 
Charbonneau (1995), namely 0.3% and 65% respectively, and 
the population size is set to 100. As it is shown, the algorithm is 
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Fig. 2. Illustration of GAlib's approach to find the global maxi- 
mum of the function (1). The different panels show the 100 indi- 
viduals of a population 0, 1,2, 3, 4, 5, 10 and 25, characterized 
by a mutation probability of 0.3% and a crossover rate of 65%. 
For the first few generations, the individuals are distributed ran- 
domly, starting from generation 5 the individuals are clearly cen- 
tered around the local maxima and continue to converge to the 
true maximum. 

capable of efficiently determining the maximum. Starting from 
generation 5 the individuals are clearly centered around the local 
maxima and continue to converge to the true maximum. At gen- 
eration 10 all individuals are very close to the global maximum, 
and we note few changes between generations 10 and 25. This 
can be explained by the large population inertia. The large popu- 
lation prohibits the fast alteration to a favorable mutation. Higher 
mutation rates could be a possible way to increase accuracy but 



we have to keep in mind that this could come at the cost of los- 
ing the global maximum (Charbonneau 1995). Furthermore it is 
not our ultimate goal to optimize this given problem in the best 
possible manner. 

The final result we get using this method after 25 genera- 
tions is: (x,y) = (0.502,0.498) with f(x,y) = 0.994. It should 
be noted that this solution is not a special case and that the algo- 
rithm delivers this result almost every time. After 1000 consec- 
utive runs we get an average result 



x = 0.501 ± 0.004 
y = 0.500 ± 0.004 



(2) 
(3) 



When we increase the number of generations, we obviously still 
recover the global maximum and reduce the standard deviation: 
for a population of 100 and 100 generations we get 



x = 0.500 ± 0.003 
y = 0.501 ± 0.003 



(4) 
(5) 



As a final example, we increase the mutation rate to 30% to show 
it improves the accuracy. Again after 100 generations, we now 
find 



x = 0.5000 ± 0.0002 
y = 0.5000 ± 0.0002 



(6) 
(7) 



Looking at Figure 2 and comparing with the result, 
f(x,y) = 0.978322 (ranked case, 40 generations), obtained by 
Charbonneau (1995), it can be seen that, for this problem, GAlib 
performs at least as good (see the 25 generation case) as the 
PIKAIA code and we can be confident to use the GAlib code 
for our purposes. 1 

2.3. FitSKIRT 

Our goal is to develop a fitting program that optimizes the pa- 
rameters of a 3D dusty galaxy model in such a way that its appar- 
ent image on the sky fits an observed image. This comes down 
to an optimization function, where the objective function to be 
minimized is the x 2 value, 



x 2 (p) 



7=1 



Imod, j(P) lobsj 



(8) 



In this expression, A^ p i x is the total number of pixels in the image, 
Imodj an d I bs j represent the flux in the j'th pixel in the simulated 
and observed image respectively, <x 7 is the uncertainty, and p 
represents the dependency on the parameters of the 3D model 
galaxy. Note that, contrary to most^ 2 problems, the uncertainty 
o-j in our case depends explicitly on the model: it can be written 
as 



obsj modj 



(P) 



(9) 



The factor cr obsj - represents the uncertainty on the flux in the j'th 
pixel of the observed image, usually dominated by a combina- 
tion of photon noise and read noise. It can be calculated using 
the so-called CCD equation (Mortara & Fowler 1981; Newberry 



1 A possible explanation for this difference might be the generational 
versus the steady state reproduction. It has to be noted, however, that 
the PIKAIA results we quote are from Charbonneau (1995). Many re- 
searchers have adapted and updated the PIKAIA code since 1995, prob- 
ably also increasing its efficiency. 
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Fig. 3. The main flowchart of the FitSKIRT procedure. Details 
on the different steps are given in the text. 



1991; Howell 2006). The factor cr mo ^j(p) represents the uncer- 
tainty on the flux in the fth pixel of the simulated image corre- 
sponding to the Monte Carlo simulation with model parameters 
p. It can be calculated during the Monte Carlo radiative transfer 
simulation according to the recipes from Gordon et al. (2001). 

The full problem we now face is to minimize the x 2 mea- 
sure (8) by choosing the best value of p in the model parameter 
space. As the link between the model image and the model pa- 
rameters is non-trivial (it involves a complete radiative transfer 
simulation), this comes down to a multidimensional, non-linear, 
non-differentiable optimization problem. 

Our approach to this problem resulted in the FitSKIRT pro- 
gram, a code that combines the radiative transfer code SKIRT 
and the genetic algorithm library GAlib. Figure 3 shows a 
flowchart of the major parts of the FitSKIRT program. 

The first step in the process consists of setting up the ingre- 
dients for the genetic algorithm. This consists in the first place 
of defining the reference image and defining the parameterized 
model describing the distribution and properties of stars and dust 
in the model galaxy. Apart from setting up this model, we also 
select the range of the parameters in the model, and define all 
genetic algorithm parameters like crossover rate, selection and 
reproduction scheme, mutation rate, etc. 

In the second phase the genetic algorithm loop (green flow) 
is started. The initial population is then created by randomly 
drawing the parameter values p from the predefined ranges so 
the genetic algorithm starts out by being uniformly spread across 
the entire parameter space. The next step is to evaluate which in- 
dividuals provide a good fit and which do not. This is done by 
starting a monochromatic SKIRT simulation using the parame- 
ters p defined by the genes of the individual. Once a simulation 
is done we compare the resulting frame and the reference im- 
age and give the corresponding individual an objective score. 
Between the creation of the image and returning the actual x 2 
value, the simulated frame is convolved with the point spread 
function (PSF) of the observed image. The third step is then con- 



cluded by returning the final x 2 (p) f° r this individual. After the 
entire population has been evaluated the best models are selected 
and offspring is created by crossover. Depending on the mutation 
rate, some of the genes of these new individuals will undergo a 
random mutation. After this step we obtain our new generation 
which is about to be evaluated next. This loop continues until a 
predefined number of generations is reached. 

Finally, when the genetic algorithm loop ends, the con- 
volved, best fitting frame is created again and the residual frame 
is determined. These residual frames are useful to investigate 
which areas of the references frames are well fitted and which 
are not. They can provide additional insight on the validity and 
consistency of the models themselves. 

In principle, the genetic algorithm searches for the best fit- 
ting model in the entire A/-dimensional parameter space, where 
N is the number of free parameters in the model. In general, 
the model image, and hence the x 2 measure (8), depends in 
a complex, non-linear way on the different parameters in the 
model, such as scalelengths, viewing angles, etc. There is one 
exception however: the radiative transfer problem is linear with 
respect to the total luminosity of the system (which always is 
one of the parameters in the model). This allows us to treat the 
total luminosity separately from the other parameters, and de- 
termine its best value outside the genetic algorithm minimiza- 
tion routine. This step effectively decreases the dimensionality 
of the parameter space the genetic algorithm has to investigate 
from N to N - 1 . In practice, this is implemented in the follow- 
ing way. Assume our 3D model is defined by the N parameters 
p = (p\, . . . Ltot). In every step of the genetic algorithm, 

the code selects a set of N - 1 parameters (pi, . . . ,Pn-i) and 
starts a SKIRT radiative transfer simulation to create a model 
image, based on these N-l parameters with a dummy value for 
L tot , and this model image is convolved with the observed image 
PSF. Before the^ 2 value corresponding to this set of parameters 
is calculated, the code determines the best fitting total luminos- 
ity of the model that minimizes the^ 2 value (8) for the particular 
values of these N-l parameters. To this aim, we can use a sim- 
ple one-dimensional optimization routine. Because no noise is 
added between the creation of the image and determining the 
final x 2 value of the genome and because the one-dimensional 
luminosity space contains no local minima, a fast and simple 
golden section search is able to determine the correct value eas- 
ily. 

Apart from finding the parameters that correspond to the 
model which best fits the observed image, we also want a mea- 
sure of the uncertainties on these parameters. A disadvantage of 
most stochastic optimization methods, or more specific genetic 
algorithms, is that such a measure is not readily available. A way 
of partially solving this problem is by using a statistical resam- 
pling technique like bootstrapping or jackknifing (Press et al. 
1992; Nesseris & Shafieloo 2010; Nesseris & Garcia-Bellido 
2012). This basically comes down to iteratively replacing a pre- 
defined number of the simulated points by the actual data points 
and comparing the resulting objective function value. Even these 
estimations of the standards errors should not be used heedlessly 
to construct confidence intervals since they remain subject to the 
structure of the data. Bootstrapping also requires the distribu- 
tion of the errors to be equal in the data set and the regression 
model (Sahinler & Topuz 2007; Press et al. 1992). Since this is 
not entirely the case for this problem, as equation (8) shows, an 
alternative approach was used to determine the uncertainty. 

Since genetic algorithms are essentially random and the pa- 
rameter space is so vast and complex, the fitting procedure can 
be repeated multiple times without resulting in the exact same 
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solution. The difficulty of differentiating between some individ- 
uals because of some closely correlated parameters will be re- 
flected in the final solutions. The entire fitting procedure used 
here consists of running five independent FitSKIRT simulations, 
all with the same optimization parameters. The standard devia- 
tion on these five solutions is set as an uncertainty on the "best" 
solution (meaning the lowest objective function value). While 
computationally expensive, this method still allows for better so- 
lutions to be found. The spread gives some insight in how well 
the fitting procedure is able to constrain some parameters and 
which parameters correlate. When the solutions are not at all co- 
herent, this can also indicate something went wrong during the 
optimization process (f.e. not enough individuals or evaluated 
generations). 

3. Application on test images 

In this section, we apply the FitSKIRT code on a simulated test 
image, in order to check the accuracy and effectiveness of this 
method in deriving the actual input parameters. 
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Fig. 4. Evolution of the determination of the dust mass in the 
FitSKIRT test simulation with one free parameter. The plot 
shows the distribution of the x 2 value for each of the 100 in- 
dividuals for 20 generations. The colors represent the different 
generations and the input value is indicated by a red line. 



3.1. The model 

The test model consists of a simple but realistic model for a dusty 
edge-on spiral galaxy, similar to the models that have been used 
to actually model observed edge-on spiral galaxies (e.g. Kylafis 
& Bahcall 1987; Xilouris et al. 1997, 1999; Bianchi 2007, 2008; 
Baes et al. 2010; Popescu et al. 2011; MacLachlan et al. 2011; 
Holwerda et al. 2012). The system consists of a stellar disc, a 
stellar bulge and a dust disc. 

The stellar disc is characterized by a double-exponential 
disc, described by the luminosity density 



J(R > z) = a 12 j ex P r tt 



exp - 



\z\ 



(10) 



with Ld the disc luminosity, h^* the radial scalelength and h z ^ 
the vertical scaleheight. For the stellar bulge, we assume the fol- 
lowing 3D distribution, 



Lb „ (m_ 
qRl U W 



where 



= Jr 2 + ^ 



(id 



(12) 



is the spheroidal radius and S n (s) is the Sersic function, de- 
fined as the normalized 3D luminosity density corresponding to 
a Sersic surface brightness profile, i.e. 



i r°° di 

n Js d/ 



(0 



dt 



I(t) = 



exp(-bt l/n ) 



(13) 



(14) 



nT(2n + 1) 

This function can only be expressed in analytical form using spe- 
cial functions (Mazure & Capelato 2002; Baes & Gentile 2011; 
Baes & van Hese 2011). Our bulge model has four free param- 
eters: the bulge luminosity Lb, the effective radius R Q , the Sersic 
index n and the intrinsic flattening q. 

The dust in the model is also distributed as a double- 
exponential disc, similar to the stellar disc, 

M d I R \ I \z\ 



Pd(R,z) = 



exp -- — 
47rh z Rd h z4 \ n RA 



exp - 



(15) 



with Md the total dust mass, and Hr A and h zA the radial scale- 
length and vertical scaleheight of the dust respectively. The cen- 
tral face-on optical depth, often used as an alternative to express 
the dust content, is easily calculated 



Cx. 



Kp(0, z) dz ■ 



kM a 
2nh l,d 



(16) 



where k is the extinction coefficient of the dust. 

In total, this 3D model has 10 free parameters, to which the 
inclination i of the galaxy with respect to the line of sight should 
be added as an 1 1th parameter. To construct our reference model, 
we selected realistic values for all parameters, based on aver- 
age properties of the stellar discs, bulges and dust discs in spiral 
galaxies (Xilouris et al. 1999; Kregel et al. 2002; Hunt et al. 
2004; Bianchi 2007; Cortese et al. 2012). The set of parameters 
is listed in the third column of Table 1 . We constructed a V-band 
image on which we will test the FitSKIRT code by running the 
SKIRT code in monochromatic mode on this model, fully tak- 
ing into account absorption and multiple anisotropic scattering. 
The optical properties of the dust were taken from Draine & Li 
(2007). The resulting input image is shown in the top panel of 
Figure 6. It is 500 x 100 pixels in size, has a pixel scale of 100 
pc/pixel and the reference frame was created using 10 7 photon 
packages. 

We now apply the FitSKIRT code with the same parameter 
model to this artificial test image and investigate whether we 
can recover the initial parameters of the model. We proceed in 
different steps, first keeping a number of parameters in the model 
fixed to their input value and gradually increasing the number of 
free parameters. 

3.2. One free parameter 

To do some first basic testing with FitSKIRT, no luminosity fit- 
ting is used on any of the images. Apart from the dust mass M d 
all parameters are fixed and set to the same of the simulated im- 
age. We let FitSKIRT search for the best fitting model with the 
dust mass ranging between 5 x 10 6 and 1.25 x 10 8 M . Note that 
this interval is not symmetric with respect to the model input 
value of 4 x 10 7 M©, which makes it slightly more difficult to 
determine the real value. 
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Table 1. Input values and fitted values of the parameters of the model used in the test simulations in Section 3. The values used to 
create the reference image (see Figure 6) are given in the third column. The fourth, fifth and sixth column list the fitted values for 
these parameters, together with their lcr error bars, for the fits with one, three and eleven free parameters respectively. 
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Fig. 5. Evolution of the dust parameters for a population of 100 individuals for the test model with three free parameters. The 
different panels show the position of each individual in the Hr^ versus h z ^ and the Hr^ versus projections of the 3D parameter 
space, for different generations (indicated by the green numbers on the top of each panel). The input values of the parameters are 
indicated by red lines. 



Figure 4 shows the evolution of a population of 100 indi- 
viduals through 20 generations with a mutation probability of 
10% and a crossover rate of 70%. In fact, when using genomes 
with only one free parameter, the crossover rate becomes quite 
meaningless. This is because the crossover between two differ- 
ent genomes will result in an offspring that is essentially exactly 
the same as the best of the parent genomes. This duplication 



will, however, still result in a faster convergence since the mu- 
tation around those genomes will generally be in a better area 
than around a random position. FitSKIRT recovers a dust mass 
of (4.00 ± 0.02) x 10 7 M Q , exactly equal to the dust mass of the 
input model. This figure shows that the parameter space is in- 
deed asymmetrical around the true value (indicated by the red 
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line) and that the entire population is gradually shifting towards 
the optimal value. 

3.3. Three free parameters 

As a second step we considere the case where we fitted three 
free parameters, more specifically the three parameters defining 
the dust distribution: the dust mass the radial scale length 
h R £ and vertical scale height h z ^. The dust parameters are hard 
to determine individually since they have to be determined di- 
rectly from the dust lane. The stellar parameters on the other 
hand are more easily determined and constrained in a larger re- 
gion outside the dust lane. The problem is close to being degen- 
erate when we look at the dust scale length and dust mass in 
exact edge-on, since changing them will roughly affect the same 
pixels. 

We consider a wide possible range for the free parameters: 
the dust scale length was searched between 3 and 9 kpc, the scale 
height between 50 and 600 pc and for the dust mass we consid- 
ered a range between 5 x 10 6 and 1.25 x 10 8 M©. In Figure 5 
we can see the evolution of a population consisting of 100 in- 
dividuals through 20 generations in the h R ^ versus h z ^ and the 
hR,d versus projections of the 3D parameter space. For now 
we only consider a crossover of 70% and disable the mutation 
so we can take a closer look at the global optimization process. 
As we can see the entire population slowly converges until, in 
the end, it is located entirely around the actual value. This is in- 
deed what we would expect from a global optimization process. 
The final best fit values are summarized in the fifth column of 
Table 1. With only three free parameters, FitSKIRT is still able 
to converge towards the exact values within the uncertainties. 
Notice that out of all parameters, the dust scale length Hrj is the 
hardest to constrain for edge-on galaxies. The values can be con- 
strained even more if we use the mutation operator because the 
search space is then investigated in a multi-dimensional normal 
distribution around these individuals. Since in the end most in- 
dividuals reside in the same area, this area will be sampled quite 
well. 

3.4. Eleven free parameters 

As a final test, FitSKIRT was used to fit all 1 1 parameters of the 
model described in section 3.1. This corresponds to a full 10+1 
dimensional optimization problem (10 parameters fitted by the 
genetic algorithm and the total luminosity determined indepen- 
dently by a golden section search). The appendix section also 
contains a comparison between genetic algorithms, Levenberg- 
Marquardt and downhill- simplex method for this specific case. 
In order to accommodate the significant increase of the param- 
eter space, we also changed the parameters of the genetic algo- 
rithm fit: we consider 100 generations of 250 individuals each 
and set the mutation rate to 5% and the crossover rate to 60%. 
We also have to keep in mind that larger populations are less 
sensitive to noise (Goldberg et al. 1991). 

The result of this fitting exercise is given in the sixth column 
of Table 1. These results are also shown graphically in Figure 6. 
The central panel shows the simulated image of our best fitting 
model, to be compared with the reference image shown on the 
top panel. The bottom panel gives the residual image. This figure 
clearly shows that the reference image is reproduced quite accu- 
rately by FitSKIRT and the residual frame shows a maximum 
of 10% discrepancy in most of the pixels, which is quite im- 
pressive considering the complexity of the problem and vastness 



of the parameter space. Apart from a good global fit, the indi- 
vidual parameters are also very well recovered: we can recover 
all fitted parameters within the Icr uncertainty interval (only the 
bulge effective radius is just outside this one standard deviation 
range). The uncertainty estimates are also useful to get more in- 
sight in which parameters are easily constrained and which are 
not. From Table 1 it is clear that the dust scale length and dust 
mass are the hardest to constrain. These parameters are close to 
degenerate when fitting edge-on spiral galaxies as the fitting rou- 
tine is most sensitive to the edge-on optical depth, i.e. the dust 
column within the plane of the galaxy. Both a dust distribution 
with a large dust mass but a small scale length and a distribution 
with a small dust mass and a large scale length can conspire to 
give similar edge-on optical depths. 

4. Application on observed data: NGC4013 

In this final section we apply FitSKIRT to recover the intrin- 
sic distribution of stars and dust in the edge-on spiral galaxy 
NGC4013. Located at a distance of about 18.6 Mpc (Willick 
et al. 1997; Russell 2002; Tully et al. 2009) and with a D 25 di- 
ameter of 5.2 arcmin, this galaxy is one of the most prominent 
edge-on spiral galaxies. It is very close to exactly edge-on and its 
dust lane can be traced to the edges of the galaxy. Conspicuous 
properties of this galaxy are its box- or peanut-shaped bulge 
(Kormendy & Illingworth 1982; Jarvis 1986; Garcia-Burillo 
et al. 1999) and the warp in the gas and stellar distribution 
(Bottema et al. 1987; Florido et al. 1991; Bottema 1995, 1996; 
Saha et al. 2009). The main reason why we selected this partic- 
ular galaxy is that it has been modelled before twice using in- 
dependent radiative transfer fitting methods, namely by Xilouris 
et al. (1999) and Bianchi (2007). This allows a direct compari- 
son of the FitSKIRT algorithm with other fitting procedures in a 
realistic context. 

We use a V-band image of NGC4013, taken with the 
Telescopio Nazionale Galileo (TNG). Full details on the obser- 
vations and data reduction can be found in Bianchi (2007). The 
image can be found in the top left panel of Figure 7. We ap- 
ply the FitSKIRT program to reproduce this V-band image, in a 
field of approximatively 5.5'xT. The same model is used as the 
test simulations discussed in Section 3, i.e. a double-exponential 
disc and a flattened Sersic bulge for the stellar distribution and 
a double-exponential disc for the dust distribution. We do not 
fix any parameters, i.e. we are again facing a 10+1 dimensional 
optimization problem. The genetic algorithm parameters are the 
same as for the test problem, i.e. 100 generations of 250 individ- 
uals each, a mutation rate of 5% and a crossover rate of 60%. 

The results of our FitSKIRT fit are given in the third col- 
umn of Table 2. We find a V-band face-on optical depth of al- 
most unity, corresponding to a dust mass of 9.9 x 10 6 M Q . In 
Figure 7 we plot the best fit to the observed V-band image of 
NGC4013 and its residuals. These two panels show that our 
FitSKIRT model provides a very satisfactory fit to the observed 
image: the residuals between the image and fit are virtually ev- 
erywhere smaller than 20%. The main exceptions where the fit is 
less accurate are the central-left region of the disc which contains 
discernible clumpy structures (which are obviously not properly 
described by our simple analytical model) and the top-right re- 
gion, which is due to the warping of the disc in NGC 4013. The 
quality of our FitSKIRT fit is quantified in Figure 8, where we 
plot the cumulative distribution of the residual values of the in- 
nermost 15,000 pixels in the fit, i.e. the number of pixels with 
a residual smaller than a certain percentage. We see that half of 
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Fig. 7. Radiative transfer model fits to a V-band image of NGC 4013. The top right panel shows the observed image, the subsequent 
rows show the best fitting SKIRT image, and the models by Bianchi (2007) and Xilouris et al. (1999). The panels on the right-hand 
side show the residual images corresponding to these models. 



the pixels have a residual less than 15%, and almost 90% of the 
pixels have a residual less than 50%. 

In the last two columns of Table 2 we list the parameters 
found by Xilouris et al. (1999) and Bianchi (2007), scaled to 
our assumed distance of 18.6 Mpc and considering the same 
value for the dust extinction coefficient k, necessary to convert 
optical depths to dust masses as in equation (16). In Figure 7 
we also compare the images and their residuals for the models 
by Xilouris et al. (1999) and Bianchi (2007) with the FitSKIRT 
model. We reconstructed these models by taking the parameters 
from Table 2 and build a V-band model image of the galaxy with 
SKIRT. The same quantitative analysis on the residual frames as 
for the FitSKIRT solution is also plotted in Figure 8. 

Looking at the model fits and their residuals, it is immedi- 
ately clear that the three models are very similar. The largest 



deviations from the observed image are found at the same lo- 
cations for the three models and hence seem to be due to im- 
possibility of our smooth analytical model to represent the true 
distribution of stars and dust in the galaxies rather than due to the 
fitting techniques. Figure 8 seems to suggest that the FitSKIRT 
solution provides a slightly better overall fit to the image, but 
this judgement might be biased by the fact that we reproduced 
the other models with our SKIRT code (their model might be 
slightly different). Looking at the model parameters in Table 2, 
it is clear that the most prominent differences between the three 
models are the different values of the face-on optical depth and 
the dust scale length. This finding is not surprising given that we 
already found in Section 3.3 and 3.4 that these parameters are 
the hardest to constrain. This also corresponds to what (Bianchi 
2007) states as the reason for the high optical depth: the large 
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Table 2. Parameters of the intrinsic distribution of stars and dust 
in NGC4013 as obtained by FitSKIRT, by Bianchi (2007) and 
by Xilouris et al. (1999). For the latter two models, the values 
are scaled to our assumed distance of 18.6 Mpc and dust masses 
are calculated using the same values for the dust extinction co- 
efficient k as used in FitSKIRT. 




r_ 0.2 0.3 

Error percentage 

Fig. 8. The cumulative distribution of the number of pixels with 
a residual smaller than a given percentage for the three model 
fits to the V-band image of NGC4013. The green, red and blue 
solid lines correspond to the FitSKIRT model, the Xilouris et al. 
(1999) model and the Bianchi (2007) model respectively. The 
dashed horizontal line corresponds to half of the total number of 
pixels. 



optical depth is compensated by the smaller scale length of the 
dust disk. 

An important aspect to take into account when comparing 
the results from the three models is the slight variations in the 
model setup and the strong differences in the radiative trans- 
fer calculations and the fitting techniques. Xilouris et al. (1999) 
used an approximate analytical method to solve the radiative 
transfer, based on the so-called method of scattered intensities 
pioneered by Henyey (1937) and subsequently perfectionized 
and applied by Kylafis & Bahcall (1987), Xilouris et al. (1997, 
1998, 1999), Baes & Dejonghe (2001a), Popescu et al. (2000, 
2011) and Mollenhoff et al. (2006). To do the fitting, they used 
a Levenberg-Marquardt algorithm. They used an approximation 
of a de Vaucouleurs model to fit the bulge and excluded the cen- 
tral 200 pc from the fit. The modelling of NGC 4013 by Bianchi 



(2007) on the other hand was based on the Monte Carlo code 
TRADING (Bianchi et al. 1996, 2000; Bianchi 2008). For the ac- 
tual fitting, he used a combination of the Levenberg-Marquardt 
algorithm (applied to model without scattering) and the amoeba 
downhill simplex algorithm. He also used an approximation of 
a de Vaucouleurs model to fit the bulge. Given these differences 
from the FitSKIRT approach, the agreement between the param- 
eters of the three models is very satisfactory. 

5. Discussion and conclusions 

In this paper, we have presented the FitSKIRT code, a tool de- 
signed to recover the spatial distribution of stars and dust from 
fitting radiative transfer models to UV, optical or near-infrared 
images. It combines a state-of-the-art Monte Carlo radiative 
transfer code with the power of genetic algorithms to perform 
the actual fitting. The noise handling properties together with 
the ability to uniformly investigate and optimize a complex pa- 
rameter space, make genetic algorithms an ideal candidate to 
use in combination with a Monte Carlo radiative transfer code. 
Using a standard but challenging optimization problem, we have 
demonstrated that the specific genetic algorithm library chosen 
for FitSKIRT, GAlib, is reliable and efficient. We could over- 
come the main shortcoming of the genetic algorithm approach, 
the lack of appropriate error bars on the model parameters, by 
deriving the spread of the individual parameters when applying 
the optimization process several times with different initial con- 
ditions. 

The FitSKIRT program was tested on an artificial reference 
image of a dusty edge-on spiral galaxy model created by the 
SKIRT radiative transfer code. This reference image was fed to 
the FitSKIRT code in a series of tests with an increasing number 
of free parameters. The reliability of the code was evaluated by 
investigating the residual frames as well as the recovery of the 
input model parameters. From both the one and three parame- 
ters fittings we concluded that the optimization process is stable 
enough and does not converge too fast towards local extrema . 
FitSKIRT recovered the input parameters very well, even for the 
full problem in which all 1 1 parameters of the input model were 
left unconstrained. 

As a final test, the FitSKIRT method was applied to deter- 
mine the physical parameters of the stellar and dust distribution 
in the edge-on spiral galaxy NGC 40 13 from a single V-band 
image. Looking at the cumulative distribution of the number of 
pixels in the residual map, we found that FitSKIRT was able to 
fit half of the pixels with a residual of less than 15% and almost 
90% of the pixels with a residual of less than 50%. Given that the 
image of NGC 4013 clearly shows a number of regions that can- 
not be reproduced by a smooth model (due to obvious clumping 
in the dust lane and a warp in the stellar and gas distribution), 
these statistics are very encouraging. The values of the fitted pa- 
rameters and the quality of the fit were compared with similar but 
completely independent radiative transfer fits done by Xilouris 
et al. (1999) and Bianchi (2007). There are some deviations be- 
tween the results, and we argue that these can be explained by 
the degeneracy between the dust scale length and face-on optical 
depth, and the differences in the model setup, radiative transfer 
technique and optimization strategy. In general, the agreement 
between the parameters of the three models is very satisfactory. 

We have demonstrated in this paper that the FitSKIRT code 
is capable of recovering the intrinsic parameters of the stellar 
and dust distribution of edge-on spiral galaxies by fitting radia- 
tive transfer models to an observed optical image. A future ex- 
tension we foresee is to include images at different wavelengths 
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in our modelling procedure. In its most obvious form, one could 
run FitSKIRT independently on different images, and use the re- 
sults obtained at different wavelengths to study the wavelength 
dependence of the stars and dust (Xilouris et al. 1997, 1998, 
1999; Bianchi 2007). In particular, the wavelength dependence 
of the optical depth, i.e. the intrinsic extinction curve, is a strong 
diagnostic for the size distribution and the composition of the 
dust grains (Greenberg & Chlewicki 1983; Desert et al. 1990; Li 
& Greenberg 1997; Weingartner & Draine 2001; Clayton et al. 
2003; Zubko et al. 2004). One step beyond this is an expansion 
of FitSKIRT in which images at different wavelengths are fit- 
ted simultaneously, and a single model is sought that provides 
the best overall fit to a set of UV/optical/NIR images. A so- 
called oligochromatic radiative transfer fitting could disentangle 
some of the degeneracies which monochromatic fitting proce- 
dures have to deal with, and can also be used to predict the cor- 
responding FIR/submm emission of the system (Popescu et al. 
2000; Baes et al. 2010; De Looze et al. 2012a,b). 

The main goal of the paper was to present the philosophy 
and ingredients behind the FitSKIRT code, and to demonstrate 
that it is capable of determining the structural parameters of 
the stellar and dust distribution in edge-on spiral galaxies using 
UV/optical/NIR imaging data. Obviously, we have the intention 
to also apply this code to real data to investigate the physical 
properties of galaxies. As the optimization process in FitSKIRT 
can cover a large parameter space and the code is almost fully au- 
tomated (which implies minimal human intervention and hence 
bias), it is a very versatile tool and is ready to be applied to a 
variety of galaxies, including larger sets of observational data. 

We are currently applying the code to a set of large edge- 
on spiral galaxies in order to recover the detailed distribution of 
stars and dust. Our main motivation is to investigate whether the 
amount of dust observed in the UV/optical/NIR window agrees 
with the dust masses derived from MIR/FIR/submm observa- 
tions, which has been subject of some debate in recent years 
(Misiriotis et al. 2001, 2006; Alton et al. 2004; Bianchi 2008; 
Baes et al. 2010; Popescu et al. 2011; MacLachlan et al. 2011). 
Major problems that have hampered much progress in this topic 
in the past were the limited wavelength coverage and sensitivity 
of the FIR/submm observations and the small number of galax- 
ies for which such observations and detailed radiative transfer 
models were available. These problems are now being cured. On 
the one hand, we now have a powerful tool available to fit op- 
tical/NIR images in an automated way without the need for and 
the bias from human intervention. On the other hand, Spitzer 
and Herschel have now provided us with deep imaging observa- 
tions of significant numbers of nearby edge-on spiral galaxies at 
various FIR/submm wavelengths (e.g. Bianchi & Xilouris 2011; 
Holwerda et al. 2012; De Looze et al. 2012a,b; Ciesla et al. 2012; 
Verstappen et al. 2012). 

We have focused our attention here on edge-on spiral galax- 
ies, which are the obvious and most popular candidates for ra- 
diative transfer modelling because of their prominent dust lanes. 
However, these are not the only targets on which FitSKIRT could 
be applied. In principle, the fitting routine can be applied on any 
geometry, but for monochromatic fitting the dust parameters can 
only be constrained for galaxies with a clear and regular sig- 
nature of dust extinction. Apart from edge-on spiral galaxies, 
an interesting class of galaxies are early-type galaxies which 
often show regular and large-scale dust lanes (e.g. Bertola & 
Galletta 1978; Hawarden et al. 1981; Ebneter & Balick 1985; 
Patil et al. 2007; Finkelman et al. 2008, 2010). We plan to ap- 
ply our FitSKIRT modelling in the future to a sample of dust- 
lane early-type galaxies, primarily focusing on those galaxies 



that have been mapped at FIR/submm wavelengths (e.g. Smith 
et al. 2012b). 

Other applications are possible and the authors welcome sug- 
gestions from interested readers. 
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Appendix A: Other optimization techniques 

As described in Section 2.2, the choice of genetic algorithms 
as our preferred optimization routine was driven by the nature 
of our problem: the minimization of a complex, numerical and 
noisy objective function in a large, multi-dimensional parame- 
ters space, characterized by several local minima. As we have 
demonstrated, the genetic algorithm approach turned out to be a 
powerful tool in this respect, and enabled us to reach our goals. 
In particular, concerning the problem of noisy objective func- 
tions, genetic algorithms are very effective: since they work on 
a population rather than iteratively on one point, that the ran- 
dom noise in the objective function is a much less determin- 
ing factor in the final solution (Metcalfe et al. 2000; Larsen & 
Humphreys 2003). On the other hand, the method is not entirely 
free of drawbacks: compared to some other optimization rou- 
tines, genetic algorithms have a relatively slow convergence rate 
on simple problems, and the error analysis is not easily treated 
in a natural way. 

To further explore other possibilities and not limit ourselves 
to this one approach only, we have performed further tests adopt- 
ing two other optimization schemes that have been used in radia- 
tive transfer fitting problems, namely the Levenberg-Marquardt 
and the downhill- simplex methods. Both methods were applied 
to the test case we have considered in Section 3.4, and their per- 
formance is compared to genetic algorithms. In the following we 
briefly describe these two routines and present the results of the 
comparison. 

A.1. Levenberg-Marquardt method 

The Levenberg-Marquardt algorithm or damped least- squares 
method (Levenberg 1944; Marquardt 1963) is an iterative op- 
timization method, that uses the gradient of the objective func- 
tion and a damping factor based on the local curvature to find 
the next step in the iteration towards the final solution. The 
Levenberg-Marquardt method is one of the most widely used 
methods for nonlinear optimization problems. It has been ap- 
plied in inverse radiative transfer modelling by Xilouris et al. 
(1999) and Bianchi (2007). The choice of this algorithm was 
possible because the radiative transfer approach used by both 
was an (approximate) analytical model (Bianchi (2007) only 
used the Levenberg-Marquardt method for fits with ray-tracing 
radiative transfer simulations where scattering was taken into ac- 
count). In that case, the gradient in the multi-dimensional param- 
eter space only requires a modest computational effort and is free 
of random noise. 

As we are dealing with a Monte Carlo radiative transfer code 
with a full inclusion of multiple anisotropic scattering, adapting 
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Fig. A.l. The function evaluations for both the Levenberg-Marquardt and downhill simplex method ranging from black (initial) to 
blue (DS) or green (LM) and the average path. The red lines indicate the values used to create the reference frame (Figure 6). Left: 
the variation of the dust mass and the scale length of the dust disk. Right: the variation of the scale length and scale height of the 
stellar disk. 




Fig. A.2. The residual frame of the reference image and the best fitting model as determined by the GA (top) and the repeatedly 
restarted Downhill Simplex method (Bottom). The parameter values are listed in Table A.l. The green line in the bottom image is 
the result of slight inaccuracies in the reproduction of the dust lane. 



this method was not trivial. In particular, the presence of Monte 
Carlo noise on the objective function evaluation makes the (nu- 
merical) calculation of the gradient difficult. In an attempt to 
overcome this problem we have combined two solutions: 

- we increase the number of photon packages in every Monte 
Carlo simulation, so that the noise is reduced to a minimum 
but leads to computationally more demanding simulations; 

- we calculate the gradient based on different different objec- 
tive function evaluations in every parameter space dimen- 
sion. This helps to mitigate both the high impact that the 
noise has on the gradient in regions close to the starting 
point, and also the unrealistic values that the gradient can 
yield when calculated on points located further away from 
the current position. 

In our implementation, we found that the number objective 
function evaluations needed to make one single Levenberg- 
Marquardt step is around 70 for our optimization problem con- 
taining 11 free parameters (Section 3.4). Keeping in mind that 
we have to use a larger number of photon packages, one sin- 
gle Levenberg-Marquardt step is comparable to approximately 
400 evaluations performed in the genetic algorithms approach, 
in terms of computational effort. 



A. 2. Downhill simplex method 

The downhill simplex method, also known as Nelder-Mead 
method or amoeba method (Nelder & Mead 1965), is another 
commonly used non-linear optimization technique. The method 
iteratively refines the search space defined by an (N + 1)- 
dimensional simplex by replacing one of its defining points, 
typically the one with the worst objective function evaluation, 
with its reflection through the centroid of the remaining N points 
defining the simplex. 

In relatively simple optimization problems, the downhill 
simplex method is known to have a relatively poor convergence, 
i.e. it is not efficient in the number of objective function evalu- 
ations necessary to find the extremum. On the positive side, the 
method only needs the evaluation of the objective function it- 
self and does not require additional information on the parame- 
ter space or the calculation of gradients. Importantly, the method 
also works better with noisy objective functions compared to 
Levenberg-Marquardt. Hence, coupling a downhill simplex rout- 
ing to a Monte Carlo radiative transfer code requires no addi- 
tional adjustments as instead needed in the case of Levenberg- 
Marquardt. The downhill simplex method was used by Bianchi 
(2007) for the general case, i.e. for fits in which scattering was 
included (and hence Monte Carlo noise was present). 
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parameter 


unit 


reference 


a 

LrA 




M d 


1U Mq 


4 


J. /j ± U.j 


A C7 i 1 ^1 

4. j / ± l.ol 


Tf 




c\ orv 


rv no i rv 1 /i 
U.98 ± U.14 


l.ZD ± U.3 / 


h R ,d 


kpc 


o.o 


COO i A TC 

5.83 ± 0.75 


^7 C\f\ i 1/1/1 

/.90 ± 1.44 




pc 


ZJl) 


24j ± 31) 


O O /I i 1 

284 ± ol 


Ltot 




1 


0.97 ± 0.09 


0.98 ± 0.13 




kpc 


4.4 


4.36 + 0.26 


4.0 ± 0.5 


K* 


pc 


500 


519 + 51 


505 ± 84 


D 1 U 




\J.JJ 




n oi + n 1 q 

U.ZJ It U. 1 V 




kpc 


2 


1.75 + 0.24 


2.16 + 0.45 


n 




2.5 


2.5 + 0.5 


3.52 ± 0.65 


q 




0.5 


0.50 + 0.04 


0.51+0.17 


i 


deg 


89 


88.9 + 0.1 


89.2 ± 0.3 



Table A.l. Input values and fitted values of the parameters of the model used in the test simulations in Section 3. The values used 
to create the reference image (see Figure 6) are given in the third column. The fourth an fifth column list the fitted values for these 
parameters, together with their Icr error bars, for the fits with eleven free parameters for the genetic algorithms and downhill simplex 
method respectively. 



A.3. Results and comparison 

When applying the Levenberg-Marquardt and downhill simplex 
algorithms to the problem discussed in Section 3.4, both meth- 
ods are initialized in the same, randomly determined position 
and new evaluations are performed until there is no significant 
improvement in the values of the objective function. The left 
panel of Figure A. 1 shows the function evaluations (points) and 
the averaged path (lines) of the Levenberg-Marquardt (green) 
and downhill simpex (blue) methods for the dust mass and 
scale length of the dust disk, two parameters which are hard to 
constrain separately. The first evaluation -or starting point- is 
marked with a black open circle, while the "real" values are in- 
dicated by the red lines. The history of the function evaluations 
is colour coded (from black at the beginning to gradually green 
or blue, for the last evaluation). 

Figure A.l visually demonstrates another fundamental com- 
plication of both the Levenberg-Marquardt and downhill simplex 
methods: they both fail to reach the area of the parameter space 
located closely to the true values. Even for parameters which 
are more easily determined, like the scale length and height of 
the stellar disk (see the right panel in Figure A.l ), the meth- 
ods do not succeed at reaching the region around the true val- 
ues. This, however, should not come as a surprise since both 
methods are known to be local rather than global optimization 
methods. This limitation can be partially overcome by repeat- 
edly restarting the algorithm in a new random initial position. 
For the Levenberg-Marquardt method, this is simply unfeasible 
as it would result in a computationally exhausting task, as the 
number of photon packages used in every Monte Carlo simu- 
lation had to be increased substantially to mitigate the effects of 
Monte Carlo noise. It was therefore not considered to be a viable 
method for this problem. 

To compare the downhill simplex method to genetic algo- 
rithms for our specific problem, we have followed the afore men- 
tioned idea, repeatedly restarting the downhill simplex optimiza- 
tion at random positions within the same boundaries, until the 
same number of function evaluations was reached. Table A.l 
lists the values of the best fit model, together with the error 
determined by the standard deviation of the five best models. 
In Figure A. 2 we show the residuals of the best fit models de- 



termined by the downhill simplex and genetic algorithms ap- 
proaches. Comparing both the residuals and the best-fit param- 
eter values, it is clear the downhill simplex method is not ca- 
pable of reproducing the nominal values with the same degree 
of accuracy reached by the genetic algorithm approach. Also, 
the spread in the parameters of the final solutions is consider- 
ably larger compared to the values obtained using the genetic 
algorithms. Again, this is the result of the "local" nature of the 
downhill simplex method. Genetic algorithms generally perform 
much better at both determining the global minimum and at han- 
dling the noise, turning out to be the most efficient and robust 
choice to solve this problem. 
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